REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data 
sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other 
aspect  of  this  collection  of  information,  including  suggestions  for  reducing  the  burden,  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information 
Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other 
provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

06/08/2014  Technical  Report  -  Final  Report 

4.  TITLE  AND  SUBTITLE 

Extended  Truncated  Hierarchical  Catmull-Clark  Subdivision 


3.  DATES  COVERED  (From  -  To) 


5a.  CONTRACT  NUMBER 


5b.  GRANT  NUMBER 


5c.  PROGRAM  ELEMENT  NUMBER 


6.  AUTHOR(S) 


5d.  PROJECT  NUMBER 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Texas  at  Austin  Austin  TX  United  States 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Office  of  Naval  Research  Arlington  VA  United  States 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 


11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

A  =  Approved  For  Public  Release  12/2/2015  No 

13.  SUPPLEMENTARY  NOTES 


14.  ABSTRACT 

In  this  paper  we  present  an  extended  Truncated  Hierarchical  Catmull-Clark  Subdivision  (eTHCCS)  method,  which 
improves  the  e  ciency  of  local  refinement  in  Truncated  Hierarchical  Catmull-Clark  Subdivision  (THCCS).  We  first 
generalize  Stam’s  Catmull-Clark  basis  functions  for  elements  with  more  than  one  extraordinary  node.  In  this  manner 
we  build  a  set  of  basis  functions  over  arbitrary  quadrilateral  meshes  and  enable  isogeometric  analysis  on  such  meshes 
without  any  preprocessing.  Then,  a  new  basis-function-insertion  scheme  is  developed  with  the  aid  of  the  truncation 

morhanicm  \A/hinh  rofinoc  nno-rinn  noinhhnrinn  olomontc  rafhor  than  txA/n-rinn  noinhhnrhnrwHc  Thorofnro  oTMCf!^ _ 

15.  SUBJECT  TERMS 


|  16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OF  RESPONSIBLE  PERSON 

a.  REPORT 

i  i 

b.  ABSTRACT 

i  i 

c.  THIS  PAGE 

i  i 

ABSTRACT 

OF 

PAGES 

U 

U 

U 

19b.  TELEPHONE  NUMBER  (Include  area  code) 

Standard  Form  298  (Rev.  8/98) 
Prescribed  by  ANSI  Std.  Z39.18 


ICES  REPORT  15-15 


June  2015 


Extended  Truncated  Hierarchical  Catmull-Clark 

Subdivision 


by 

Xiaodong  Wei,  Yongjie  Jessica  Zhang,  Thomas  J.R.  Hughes,  and  Michael  A.  Scott 


The  Institute  for  Computational  Engineering  and  Sciences 

The  University  of  Texas  at  Austin 
Austin,  Texas  78712 


Reference:  Xiaodong  Wei,  Yongjie  Jessica  Zhang,  Thomas  J.R.  Hughes,  and  Michael  A.  Scott,  "Extended 
Truncated  Hierarchical  Catmull-Clark  Subdivision, "  ICES  REPORT  15-15,  The  Institute  for  Computational 
Engineering  and  Sciences,  The  University  of  Texas  at  Austin,  June  2015. 


Extended  Truncated  Hierarchical  Catmull-Clark  Subdivision 


Xiaodong  Wei3,  Yongjie  Jessica  Zhang3’*,  Thomas  J.R.  Hughesb,  Michael  A.  Scottc 

aDepartment  of  Mechanical  Engineering,  Carnegie  Mellon  University,  Pittsburgh,  PA  15213,  USA 
b Institute  for  Computational  Engineering  and  Sciences,  The  University  of  Texas  at  Austin,  Austin,  TX  78712,  USA 
c Department  of  Civil  and  Environmental  Engineering,  Brigham  Young  University,  Provo,  UT  84602,  USA 


Abstract 

In  this  paper  we  present  an  extended  Truncated  Hierarchical  Catmull-Clark  Subdivision  (eTHCCS)  method,  which 
improves  the  efficiency  of  local  refinement  in  Truncated  Hierarchical  Catmull-Clark  Subdivision  (THCCS).  We  first 
generalize  Stam’s  Catmull-Clark  basis  functions  for  elements  with  more  than  one  extraordinary  node.  In  this  manner 
we  build  a  set  of  basis  functions  over  arbitrary  quadrilateral  meshes  and  enable  isogeometric  analysis  on  such  meshes 
without  any  preprocessing.  Then,  a  new  basis-function-insertion  scheme  is  developed  with  the  aid  of  the  truncation 
mechanism,  which  refines  one-ring  neighboring  elements  rather  than  two-ring  neighborhoods.  Therefore,  eTHCCS 
significantly  improves  the  efficiency  of  local  refinement  compared  with  THCCS,  as  demonstrated  by  one  benchmark 
problem  and  several  complex  models.  Moreover,  eTHCCS  is  also  proved  to  preserve  the  input  geometry  and  produce 
nested  spaces. 

Keywords :  Local  refinement,  truncated  hierarchical  Catmull-Clark  subdivision,  arbitrary  control  meshes, 
isogeoemtric  analysis. 


1.  Introduction 

Local  refinement  and  arbitrary  topology  have  been  two  key  concerns  since  isogeometric  analysis  [12,  5]  was  pro¬ 
posed.  Isogeometric  analysis  aims  to  integrate  engineering  design  with  simulation  by  employing  the  basis  used  for  ge¬ 
ometric  design  in  analysis.  Local  refinement  was  studied  in  hierarchical  B-splines  [8,  14,  25,  2],  T-splines  [23,  22,  1], 
PHT-splines  [6],  LR-splines  [7,  13],  THB-splines  [9,  10,  29]  and  Truncated  Hierarchical  Catmull-Clark  Subdivision 
(THCCS)  [28].  Among  them,  hierarchical  B-splines  and  THB-splines  are  restricted  to  a  topologically  rectangu¬ 
lar  domain,  and  thus  they  have  limited  applications  in  complex  geometries.  T-splines  support  arbitrary  topologies 
[27,  26,  16],  but  the  original  development  of  T-splines  may  violate  one  prerequisite  requirement  in  analysis,  i.e.,  linear 
independence  [3].  A  mildly  restricted  subset,  analysis- suitable  T-splines  [15,  21],  was  subsequently  developed  to  ad¬ 
dress  the  linear  independence  issue.  However,  local  refinement  of  analysis-suitable  T-splines  may  propagate  beyond 
the  domain  of  interest  via  excessive  T-junction  extension.  PHT-splines  support  arbitrary  topology  and  element-wise 
local  refinement  with  a  trade-off  in  almost  twice  the  number  of  degrees  of  freedom  and  reduced  continuity  (C1  for 
cubic  splines).  THCCS,  on  the  other  hand,  addresses  both  local  refinement  and  arbitrary  topology  without  introducing 
extra  degrees  of  freedom  and  attains  C2  continuity. 

THCCS  [28]  was  developed  based  on  Catmull-Clark  subdivision  [4,  30,  19].  Catmull-Clark  subdivision  is  a 
popular  quadrilateral-based  subdivision  scheme  that  is  generalized  from  mid-knot  insertion  of  bi-cubic  B-splines 
to  arbitrary  topologies.  A  Catmull-Clark  surface  is  represented  by  a  quadrilateral  mesh,  which  is  obtained  by  an 
iterative  and  global  refinement  of  an  initial  coarse  quadrilateral  mesh.  The  sequence  of  refined  meshes  converges  to 
a  limit  surface  that  is  C2 -continuous  everywhere  except  C1  -continuous  at  extraordinary  vertices  (a  vertex  is  called 
an  extraordinary  vertex  if  it  has  other  than  four  quadrilaterals  adjacent  to  it,  and  it  is  a  regular  vertex  otherwise).  In 
this  paper,  vertex  and  quadrilateral  are  used  in  the  physical  domain,  and  correspondingly  in  the  parametric  domain, 
their  preimages  are  referred  to  as  node  and  element ,  respectively.  Alternatively,  the  limit  surface  can  be  directly 
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evaluated  from  the  initial  quadrilateral  mesh  by  Stam’s  basis  functions  [24].  The  quadrilateral  mesh  is  called  the 
control  mesh  of  the  Catmull-Clark  surface.  Based  on  Stam’s  Catmull-Clark  basis  functions,  THCCS  takes  advantage 
of  the  hierarchical  structure  and  truncation  mechanism  to  support  local  refinement  and  preserve  geometry.  THCCS  has 
desired  properties  for  geometric  design  and  analysis  such  as  partition  of  unity,  convex  hull,  linear  independence,  and 
support  of  arbitrary  topology.  However,  two  requirements  in  THCCS  restrict  the  efficiency  of  local  refinement:  (1)  at 
most  one  extraordinary  vertex  is  allowed  in  each  quadrilateral  in  the  control  mesh;  and  (2)  for  each  to-be-refined  basis 
function,  all  its  two-ring  neighboring  elements  have  to  be  refined.  Requirement  (1)  is  inherited  from  the  development 
of  Stam’s  Catmull-Clark  basis  functions  [24].  In  THCCS,  this  requirement  is  satisfied  by  refining  all  the  quadrilaterals 
that  contain  more  than  one  extraordinary  control  point  (such  quadrilaterals  are  called  invalid).  The  input  quadrilateral 
mesh  may  contain  a  large  number  of  quadrilaterals  that  violate  Requirement  (1),  leading  to  almost  global  refinement. 
Requirement  (2)  follows  the  basis-function-refinement  manner  in  hierarchical  B-splines  [14,  11,  25,  2].  For  cubic 
hierarchical  B-splines  and  Catmull-Clark  subdivision,  however,  such  basis-function-refinement  needs  to  refine  all  the 
two-ring  neighboring  elements  of  each  to-be-refined  basis  function,  which  is  not  efficient  in  capturing  abrupt  changes 
in  geometric  or  solution  features. 

To  improve  the  efficiency  of  local  refinement  in  THCCS,  in  this  paper  we  develop  the  extended  Truncated  Hierar¬ 
chical  Catmull-Clark  Subdivision  (eTHCCS).  We  first  make  a  straightforward  generalization  of  Stam’s  Catmull-Clark 
basis  functions  to  enable  their  direct  application  on  arbitrary  control  meshes,  where  more  than  one  vertex  is  allowed  in 
a  single  quadrilateral.  This  generalization  eliminates  the  requirement  of  refining  invalid  quadrilaterals,  and  provides  a 
set  of  basis  functions  for  isogeometric  analysis  on  arbitrary  quadrilateral  meshes.  The  analysis-suitability  on  irregular 
quadrilateral  meshes  is  also  studied  for  the  generalized  Catmull-Clark  basis  functions.  The  main  contribution  of  this 
paper  is  the  development  of  a  new  basis -function-insertion  scheme  for  local  refinement  to  improve  its  efficiency  by 
releasing  the  restriction  on  the  to-be-refined  region.  In  THB -splines  or  THCCS,  the  to-be-refined  region  is  restricted  to 
be  the  support  of  to-be-refined  basis  functions,  which  is  the  union  of  two-ring  neighboring  elements  of  the  associated 
nodes  in  the  case  of  cubic  splines.  In  addition  to  preserving  all  the  properties  of  THCCS,  eTHCCS  therefore  achieves 
a  higher  efficiency  in  local  refinement. 

The  remainder  of  this  paper  is  organized  as  follows.  Section  2  reviews  the  basic  concepts  of  THCCS.  Section 
3  introduces  the  detailed  development  of  eTHCCS.  In  Section  4,  we  study  the  geometry  preservation  and  nested 
property  of  eTHCCS.  We  then  demonstrate  the  efficiency  of  the  proposed  method  in  Section  6,  and  conclude  the 
paper  in  Section  7. 

2.  A  Review  of  Truncated  Hierarchical  Catmull-Clark  Subdivision 

In  this  section  we  briefly  review  the  key  concepts  of  THCCS.  For  details  one  may  refer  to  related  literature 
[4,30,24,14,25,2,9,28]. 

2.7.  Stam ’s  Catmull-Clark  Basis  Functions 

THCCS  was  developed  based  on  Stam’s  Catmull-Clark  basis  functions.  It  is  a  natural  starting  point  to  introduce 
Catmull-Clark  subdivision  and  its  basis  functions.  Catmull-Clark  is  one  of  the  most  popular  subdivision  schemes 
used  in  the  CAD  community.  The  refinement  (or  subdivision)  in  Catmull-Clark  subdivision  is  generalized  from  mid¬ 
knot  insertion  of  bi-cubic  B-splines  [4].  Each  new  vertex  in  the  refined  mesh  is  calculated  via  a  weighted  average 
of  neighboring  vertices  in  the  original  mesh.  This  linear  relationship  can  be  expressed  by  a  so-called  subdivision 
matrix.  The  repeated  global  refinement  generates  a  sequence  of  meshes,  Af°, . . . ,  Af\  where  Af°  is  the  initial  input 
quadrilateral  mesh,  and  n  is  the  number  of  subdivisions.  As  n  goes  to  infinity,  Mn  converges  to  a  limit  surface.  An 
alternative  way  to  obtain  the  limit  surface  takes  advantage  of  the  Stam’s  basis  functions  [24].  These  basis  functions 
are  analogous  to  B -spline  basis  functions,  whereas  each  mesh  M?  (0  <  b  <  n)  behaves  as  a  control  mesh.  Thus  we 
can  express  the  limit  surface  <Snmit  by  a  mapping  from  the  parametric  domain  to  the  physical  domain, 

Ne 
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where  #^(£,77)  are  Stam’s  Catmull-Clark  basis  functions,  £  and  77  are  parametric  coordinates,  Pj  are  control  points 
(or  vertices)  in  the  physical  domain,  and  Ne  is  the  number  of  basis  functions  of  M* .  We  call  £  the  subdivision  level. 
Without  loss  of  generality,  we  introduce  the  basis  functions  at  Level  £  as  follows. 
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Figure  1:  Local  parametric  elements,  (a)  A  regular  element  (blue);  (b)  an  irregular  element  (blue);  (c)  subdivision  of  (b);  and  (d)  subdivisions  in 
the  parametric  domain.  The  local  parametric  domain  is  [0, 1]  x  [0, 1].  Nodes  are  locally  numbered  with  respect  to  the  element  marked  in  blue. 


A  Catmull-Clark  subdivision  surface  generally  does  not  have  a  global  rectangular  parametric  domain  due  to  the 
presence  of  extraordinary  nodes.  A  parametric  element  is  locally  associated  with  each  quadrilateral  in  the  control  mesh 
and  each  control  point  has  a  corresponding  node  in  the  parametric  domain.  If  an  element  contains  any  extraordinary 
node,  it  is  irregular,  and  otherwise  it  is  regular.  Fig.  1(a)  shows  a  local  mesh  surrounding  a  regular  element  (shaded  in 
blue).  A  total  of  16  basis  functions  have  support  over  this  regular  element  because  a  Catmull-Clark  basis  function  has 
local  support  over  its  two-ring  neighboring  elements.  They  are  actually  bi-cubic  uniform  B -spline  basis  functions, 

n)  =  i  =  1, 2, ....  16,  (2) 

where  and  “/”  represent  the  remainder  and  division,  respectively,  and  for  t  e  [0, 1]  we  have 

b0(t)  =  (1  -  3 1  +  3 12  -  t3)/ 6,  bi{t)  =  (4  -  6t2  +  3f3)/6,  (3) 

b2(t)  =  (1  +  3 1  +  3 12  -  3 f3)/6,  b3(t)  =  t3 / 6.  (4) 

Fig.  1(b)  shows  a  local  mesh  surrounding  a  valence-3  vertex1,  where  an  irregular  element  (Qq)  is  marked  in 
blue.  The  surrounding  2 N  +  8  (N  is  the  valence  number  and  here  N  =  3)  basis  functions  have  support  on  Qq,  whose 
associated  vertices  are  locally  labeled  in  the  manner  shown  in  Fig.  1(b).  The  2 N  +  8  basis  functions  over  are 
derived  by  infinitely  subdividing  [24] .  In  the  first  subdivision,  is  subdivided  into  one  smaller  irregular  element 
Qq+1  and  three  regular  elements  Qfk+l  (k  =  1, 2,  3);  see  Fig.  1(c).  The  new  2N  + 17  vertices  in  Fig.  1(c)  can  be  obtained 
from  the  IN  +  8  vertices  in  Fig.  1(b)  via  a  subdivision  matrix  A(2^v+i7)x(2v+8)-  As  the  basis  functions  are  well-defined 
on  Q?k+1  (k  =  1, 2,  3)  as  in  Fig.  1(a),  the  limit  surface  corresponding  to  the  3/4  parametric  domain  of  Df0  is  represented 
as 


(Bf)Tpf  =  (Bc+y  P£+i  =  (Bc+i)' AP"  =  (A1  BC+LY  P€, 
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where  =  [B{, . . . ,  ^+8]r,  P €  =  Pe2N+ 8\T,  Bm  =  [B{+\  ...,Bl 

Eq.  (5)  holds  for  any  P6,  so  we  have 
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2N+\1 
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B^  =  ArB^+1.  (6) 

The  remaining  1/4  parametric  domain  of  is  the  irregular  element  ^q+1*  nee(^  t0  father  subdivide  Q^+1  to 
find  another  3/4  parametric  domain  of  where  basis  functions  are  well-defined.  By  repeating  this  procedure,  the 


!The  valence  number  of  a  node  is  the  number  of  elements  adjacent  to  it. 
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parametric  domain  corresponding  to  the  irregular  element  becomes  exponentially  smaller,  as  shown  in  Fig.  1(d).  The 
repeated  subdivision  occurs  in  the  irregular  element.  Note  that  the  subdivision  of  Q^+1  only  involves  the  first  2N  +  8 
vertices  in  Fig.  1(c).  Therefore  the  sub-matrix  consisting  of  the  first  2N  +  8  rows  of  A  will  be  repeatedly  used,  denoted 
as  A(2a/’+8)x(2v+8)  •  For  computational  efficiency,  the  eigenstructure  (A,  V)  of  A  (AV  =  VA)  is  employed  such  that  only 
diagonal  matrix  multiplication  is  required.  In  this  manner,  the  basis  functions  at  Level  £  over  the  irregular  element 
is  derived  as 

n)  =  (V_l),A"_l(P/-AV)/b(4r,  if),  (7) 

where  b(£,  rf)  are  the  uniform  bi-cubic  B -spline  basis  functions  as  in  Eq.  (2).  Given  parametric  coordinates  (£,  rf),  we 
perform  subdivision  n  times  to  restrict  (£,  rj)  into  a  regular  element  (Q",  k  -  1, 2, 3)  as  in  Fig.  1(d).  P k  (k  =  1,2, 3) 
is  a  selection  matrix  to  locate  such  regular  element.  The  configuration  of  A  and  A  can  be  found  in  the  appendix  A  in 
S tarn’s  work  [24],  which  only  depends  on  the  valence  of  the  extraordinary  node.  When  the  parametric  values  (£,  rj) 
approach  zero,  Eq.  (7)  is  defined  as  a  limit  case  where  An~l  becomes  a  matrix  such  that  only  its  first  diagonal  element 
is  non-zero.  This  is  because  in  the  diagonal  matrix  A,  all  the  elements  are  positive  and  smaller  than  1  except  for  the 
first  element,  which  equals  to  1.  Eq.  (6)  is  also  valid  over  Q^+1  [28],  where  basis  functions  are  defined  by  Eq.  (7). 
Eq.  (6)  indicates  a  general  relationship  between  basis  functions  at  two  consecutive  levels.  We  call  this  relationship 
refinability ,  and  we  define  high-level  basis  functions  B^+1  as  the  children  of  low-level  basis  functions  in  Eq.  (6). 
Refinability  is  fundamental  to  the  construction  of  THCCS. 

However,  Eq.  (7)  does  not  work  for  all  quadrilateral  meshes  because  its  derivation  requires  that  each  quadrilateral 
contains  at  most  one  extraordinary  node.  Usually,  an  input  quadrilateral  mesh  needs  to  be  globally  refined  once  before 
S tarn’s  basis  functions  are  applied. 

2.2.  Truncated  Hierarchical  Catmull- Clark  Subdivision 

THCCS  [28]  generalizes  truncated  hierarchical  B-splines  (THB-splines)  [9]  to  control  meshes  of  arbitrary  topol¬ 
ogy.  THB- splines  were  developed  to  further  modify  the  hierarchical  B- spline  basis  functions  to  form  a  partition  of 
unity  and  to  decrease  overlapping.  However,  THB-splines  can  only  be  used  to  represent  the  geometries  with  a  global 
parametric  domain.  Complex  geometries  unavoidably  involve  extraordinary  nodes  and  thus  cannot  be  mapped  onto 
a  global  parametric  domain.  Based  on  this  fact,  THCCS  was  developed  as  an  attempt  to  address  local  refinement  on 
arbitrary  topologies. 

Starting  from  a  valid  input  quadrilateral  mesh,  we  recursively  construct  THCCS.  The  recursive  manner  allows  us 
to  consider  two  consecutive  levels  at  one  time.  We  now  construct  Level  £  +  1  from  Level  £.  Levels  basis  functions 
are  denoted  as  & .  The  THCCS  space  can  be  enlarged  by  replacing  the  identified  basis  functions  (0r  c  Se)  with  their 
children  (chd &r).  After  the  refinement  of  &r,  we  define  active  Levels  basis  functions  as  and  the  children 

basis  functions  of  B,  (chd^)  are  the  active  Level-(^  +1)  basis  functions.  Only  active  basis  functions  are  collected 
into  THCCS  basis.  Besides,  if  a  Level-^  basis  function  B f  €  Bl\B[  has  any  children  contained  in  chdiSf,  it  has 
to  be  truncated  in  order  to  form  a  partition  of  unity  and  to  preserve  the  geometry.  The  truncation  is  performed  by 
discarding  active  children  from  the  refinability  relationship.  Recall  that  according  to  refinability,  Bj  can  be  expressed 
as  Bl{  -  YjbmecMb{  cij^j+1  ’  w^ere  coefficients  ctj  come  from  the  subdivision  matrix  A.  Then  the  truncation  of  B\  is 
obtained  by 


tmnflf  =  L  (8) 

ZU1  echdif  and  Z^+1g chdSf 

Finally,  we  collect  active  basis  functions  at  Level  £  and  Level  £  +  1  into  the  THCCS  basis.  Other  high  levels  can  be 
constructed  likewise  and  the  THCCS  basis  is  updated  accordingly.  The  above  procedure  can  be  recursively  performed 
until  the  desired  maximum  level  £mSiX  is  reached. 

THCCS  basis  functions  satisfy  partition  of  unity,  convex  hull  property,  linear  independence  with  the  support  of 
local  refinement,  and  arbitrary  topologies.  However,  the  basis-function-refinement  nature  requires  refinement  of  two- 
ring  neighborhood  of  a  basis  function  at  each  refinement  step,  leading  to  the  refinement  of  all  the  elements  within 
its  support.  eTHCCS  is  therefore  developed  to  improve  the  efficiency  of  local  refinement  with  a  new  basis-function- 
insertion,  in  which  we  only  need  to  refine  one-ring  neighboring  elements.  This  new  scheme  releases  the  strong 
restriction  of  the  minimum  number  of  to-be-refined  elements.  Details  of  eTHCCS  will  be  discussed  next  in  Section  3. 
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3.  Development  of  eTHCCS 


In  this  section,  we  discuss  how  to  develop  eTHCCS.  First,  we  generalize  the  Catmull-Clark  basis  functions  for 
elements  with  arbitrary  numbers  of  extraordinary  nodes,  eliminating  the  requirement  of  refining  such  elements  in 
order  to  use  S tarn’s  basis  functions.  In  this  manner  we  build  the  basis  suitable  for  isogeometric  analysis  over  arbitrary 
quadrilateral  meshes.  Then  we  develop  a  new  basis-function-insertion  scheme  to  improve  the  locality  of  refinement, 
releasing  the  restriction  of  the  to-be-refined  region.  With  this  scheme,  we  can  refine  even  only  one  element  at  each 
refinement  step,  rather  than  refining  two-ring  neighborhoods  of  elements  in  THCCS.  Finally,  we  use  the  generalized 
Catmull-Clark  basis  functions  and  basis-function-insertion  scheme  to  construct  eTHCCS. 

3.1.  Generalized  Catmull-Clark  Basis  Functions 

Here  we  aim  to  generalize  the  derivation  of  Catmull-Clark  basis  functions  to  invalid  elements.  Recall  that  an 
invalid  element  contains  more  than  one  extraordinary  node.  Once  an  invalid  element  is  subdivided,  the  resulting 
four  high-level  elements  are  all  valid  over  which  basis  functions  are  defined  in  either  Eq.  (2)  or  Eq.  (7).  Instead  of 
subdividing  invalid  elements  first  and  then  applying  the  subdivision  matrix  A  to  derive  Eq.  (7),  in  the  following  we 
introduce  a  more  general  subdivision  matrix,  denoted  as  S,  and  directly  apply  it  to  derive  generalized  Catmull-Clark 
basis  functions  over  an  invalid  element. 


(a)  (b) 


Figure  2:  Generalized  subdivision,  (a)  An  invalid  element  Q®  with  its  two-ring  neighboring  nodes  labeled;  and  (b)  subdividing  Q®  into  four 
high-level  elements  (k  =  0, . . . ,  3),  whose  two-ring  neighboring  nodes  are  labeled  in  the  similar  manner. 

Given  an  invalid  element,  let  us  first  locally  label  its  two-ring  neighboring  nodes.  We  follow  the  manner  of  labeling 
as  in  Fig.  1(b),  labeling  the  one-ring  neighboring  nodes  of  the  extraordinary  node  clockwise.  Note  that  any  corner 
node  of  an  invalid  element  can  be  an  extraordinary  node.  For  instance,  in  Fig.  2(a),  is  the  invalid  element  under 
study.  Let  No,  N\,  N2  and  A3  be  the  valence  numbers  of  four  comer  nodes  of  Q®,  labeled  as  Node  1,  Node  6,  Node  5 
and  Node  4,  respectively.  We  start  labeling  the  one-ring  neighboring  nodes  of  Node  1  in  clockwise,  from  2,  3,  until 
i0  =  2Ao  +  1.  Then  we  consider  the  one-ring  neighboring  nodes  of  Node  6.  Among  them  only  2N\  +  1  -  6  =  2N\  -  5 
nodes  remain  unlabeled,  and  we  label  them  from  z'o  +  1  until  i\ ,  where  i\  =  z'o  +  2Afi  -  5  =  2(Ao  +  N\)  -  4.  Here  Node 
6  is  valence-3,  so  we  have  N\  =  3  and  i\  =  z0  +  1  •  Next  for  Node  5,  there  are  2N2  +  1  -  6  =  2N2  -  5  nodes  remaining 
unlabeled  if  A3  >  3;  otherwise,  there  are  2N2  +  1  -  7  =  2N2  -  6  unlabeled  nodes  if  A3  =  3,  because  one  more  node  was 
labeled  already  when  we  labeled  the  one-ring  neighboring  nodes  for  Node  1.  We  label  them  from  i\  +  1  until  h  where 
i2  =  h  +  2A2  -  5  =  2(A0  +  Ni  +  N2)  -  9  (N3  >  3)  or  i2  =  h  +2  N2-6  =  2  (A0  +  Nx  +  N2)  -  10  (A3  =  3).  Similarly  for 
Node  4,  the  number  of  unlabeled  nodes  is  2N3  +  1  -  7  =  2N3  -  6  if  A3  =  3,  otherwise  it  is  2N3  +  1  -  8  =  2N3  -  7.  We 
label  them  from  i2  +  1  until  i3  where  i3  =  2  A)  -  16.  The  labels  are  shown  in  detail  in  Fig.  2(a).  The  associated 
basis  functions  are  the  generalized  basis  functions  to  be  derived,  and  we  denote  them  in  the  vector  form, 

B  =  [BuB2,...,Bh]T.  (9) 
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Correspondingly,  their  control  points  are  denotes  as  P  =  [Pi,  P2, . . . ,  P*3]r.  As  shown  in  Fig.  2(b),  we  subdivide 
Dq  into  four  smaller  elements  (k  =  0, . . . ,  3)  and  follow  in  the  same  manner  to  label  their  two-ring  neighboring 
nodes.  Note  that  Node  1,  Node  jo  +  1 ,  Node  j\  +  1  and  Node  72  +  1  have  the  same  valence  numbers  as  the  four  corner 
nodes  of  in  Fig.  2(a).  Similarly,  we  can  derive  the  following  relationships,  jo  =  2No  +  1,  j\  =  2{No  +  N\)  -  1, 
72  =  2{No+N\  -\-N2)-2  and  73  =  2(Ao+Afi  -fA^+A^)^.  The  corresponding  basis  functions  are  Stam’s  Catmull-Clark 
basis  functions  defined  by  either  Eq.  (2)  or  Eq.  (7),  denoted  as 

B  =  [BuB2,...,Bh]T.  (10) 

Their  control  points  are  P  =  [P\,P2,  ■  ■  ■ ,  PJ3  \T ,  which  are  calculated  by 

P  =  SP,  (11) 

where  S  can  be  directly  obtained  from  the  Catmull-Clark  subdivision  rule.  For  instance,  Vertex  1  is  relocated  to  a 
weighted  average  of  its  neighboring  vertices,  whose  indices  are  I  =  {1, 2, 3, . . . ,  i'o}.  Then  the  1st  row,  /- th  (  j  e  7) 
column  element  of  S  is  filled  with  the  corresponding  Catmull-Clark  subdivision  coefficient.  Other  elements  in  S  can 
be  filled  in  the  same  manner.  The  configuration  of  S  is  given  in  Appendix  A.  Assume  that  the  evaluation  of  Q®  using 
B  yields  the  same  limit  surface  as  that  of  Q,1,  ~  Qi  using  B,  and  we  have 

BrP  =  BrP  =  BrSP  =  (SrB)rP.  (12) 

Eq.  (12)  holds  no  matter  what  the  values  of  P  are.  Therefore  we  have 

B  =  SrB.  (13) 


j2+1 


4_j2+1_j2+2 


Jo+ 

h Z 


] k 


3  jif4 


dl?Z2fcjl: 


jol3  h+2 


+1  jif5 
h 


(c)  k  =  2 


Figure  3:  Generalized  selections  of  basis  functions  for  high-level  elements  Q®  (k  =  0, . . . ,  3). 


We  next  find  the  relationship  between  B  and  Eq.  (2)  or  Eq.  (7)  so  that  we  can  obtain  an  explicit  expression 
for  B.  Given  a  pair  of  parametric  coordinates  (£,  rj)  in  we  can  locate  it  in  one  of  the  four  elements  (Q^,  where 
k  =  0, . . . ,  3)  in  Fig.  2(b),  and  only  the  two-ring  basis  functions  of  that  element  are  non-zero  at  (£,  rj).  Thus  B  is 
defined  piecewise.  For  instance,  if  0  <  £  <  1/2  and  0  <  77  <  1/2,  (£,  rj )  is  located  in  the  irregular  element  and  only 
the  two-ring  basis  functions  are  non-zero  over  as  shown  in  Fig.  3(a).  These  2No  +  8  basis  functions  (denoted  as 
Bo)  are  selected  from  B  in  Eq.  (13).  To  directly  use  Eq.  (7),  we  need  to  sort  Bo  such  that  the  basis  functions  in  Bo 
have  the  same  order  as  in  Fig.  1(b),  that  is, 

Bo  =  [Bi,B2,  .  .  .  ,P7o,P7l  +  i,P7o+2,P7o+l,P21,P21+3,P22  +  l,Pj2+3]r-  (14) 

Note  that  in  Bo,  the  first  2No  +  1  basis  functions  correspond  to  the  extraordinary  node  and  its  one-ring  neighboring 
nodes  (from  1  to  jo),  and  the  next  4  basis  functions  are  sorted  along  the  opposite  rj 0  direction,  and  the  last  3  basis 
functions  follow  the  opposite  £0  direction.  This  specific  manner  of  sorting  follows  Stam’s  [24].  Then  we  can  directly 
obtain  Bo  using  Eq.  (7)  by  replacing  B^  ( l  -  0)  with  Bo.  We  define  a  set  of  pairs  as 

Vo  ={(1, 1),  (2, 2), . . . ,  (2A0  +  1, 70),  (2No  +  2,  j,  +  1),  (2 N0  +  3, 70  +  2),  (2 N0  +  4, 70  +  1),  (2Ao  +  5, 71), 

(2No  +  6, 71  +  3),  (2No  +  7, 72  +  1),  (2  N0  +  8, 72  +  3)}, 
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which  represents  the  element  correspondence  between  B0  and  B.  For  (/,  j)  €  Po ,  the  /- th  element  in  B0  is  the  7-th 
element  in  B.  For  instance  in  (2,/Vo  +  2,  j\  +  1),  the  (2Ao  +  2)-th  element  in  Bo  is  Bjl+ 1,  which  is  also  the  (j\  +  l)-th 
element  in  B.  Moreover  since  0  <  £  <  1/2  and  0  <  77  <  1/2,  the  basis  functions  in  B  other  than  those  in  Bo  are  all  zero. 
Then  we  can  obtain  B  by  a  linear  transformation,  B  =  PqB0,  where  Po  is  the  so-called  selection  matrix  with  respect 
to  Dq  and  it  maps  the  selected  basis  functions  to  the  correct  positions  in  B.  Po  is  a  matrix  of  dimension  (2A/o  +  8)  x  73 
and  all  its  elements  are  zero  except  that  the  /-th  row,  7-th  column  element  is  1,  where  (/,  7)  e  Po.  Plugging  B  =  P^Bo 
into  Eq.  (13),  we  have 

6  =  (P0S)rB0,  (16) 

when  0  <  £  <  1/2  and  0  <  77  <  1/2. 

Likewise,  other  selections  are  shown  in  Fig.  3(b-d)  for  (£,77)  e  [1/2, 1]  x  [0, 1/2),  (£,77)  e  [1/2, 1]  x  [1/2, 1] 
and  (£,  77)  e  [0, 1/2)  x  [1/2, 1],  respectively.  The  corresponding  selection  matrix  P^  has  the  dimension  (2A^  +  8)  x  73 
(k  =  1, 2, 3)  and  the  selected  basis  functions  B^  are  sorted  with  the  aid  of  the  corresponding  set  of  pairs  P £.  Note  that 
Q}  in  Fig.  3(b)  is  also  an  irregular  element,  so  Pi  and  Bi  are  similar  to  Po  and  Bo,  respectively.  However  in  Fig.  3(c, 
d),  Q}2  and  D*  are  regular  elements,  and  thus  the  selected  basis  functions  B2  and  B3  can  be  directly  obtained  from  Eq. 
(2).  B2  and  B3  are  sorted  following  Fig.  1(a).  For  example,  the  set  of  pairs  of  B2  is  defined  as 

r2  ={(1,  h\  (2, 71  +  5),  (3,  71  +  4),  (4, 72  +  2),  (5, 71  +  2),  (6, 71  +  1),  (7, 71  +  3),  (8, 72  +  1), 

(9, 70  +  3),  (10, 70  +  2),  (11,5),  (12, 4),  (13, 7 0,  (14, 70  +  D,  (15, 6),  (16, 1)}.  (  ' 1 

In  summary,  we  derive  the  basis  functions  over  an  invalid  element  as 


8(6  TJ)  =  (P*S)rBt(&,  m),  k  =  0, 1, 2, 3, 


OB) 


where  (6,  Vi<)  are  the  parametric  values  defined  in  the  local  coordinate  system  of  £2]  (k  =  0, . . . ,  3).  Note  that  in  Fig. 
3,  the  local  coordinate  systems  of  four  elements  Qj.  =  [0, 1]  x  [0, 1]  (k  -  0, . . . ,  3)  do  not  coincide  with  that  of  D[j. 
Given  the  parametric  coordinates  (6  >/)  in  Q®,  we  need  to  transform  them  to  parametric  values  consistent  with  the 
local  coordinate  system  of  flj. ,  where  k  =  0, . . . ,  3.  We  have 


:  (£o,  no ) 

n\  ■  (turn) 
n'2  :  (&,m) 
fl'  :  i£o,no) 


=  (26, 2/7) 

=  (277, 2  -2^) 

=  (2  -  2£  2  -  2rf) 

=  (2-2i7,26 


if  0  <  4  <  \ 
if  ±<£<1 
if  ^<^<1 
if  0  <  f  <  \ 


and  0  <  77  < 
and  0  <  77  <  2, 
and  \  <  77  <  1, 
and  4  <rj  <  1. 


(19) 


Remark  3.1.  In  classification,  we  have  a  total  of  three  types  of  elements  in  the  control  mesh:  regular  elements, 
irregular  elements  and  invalid  elements.  The  first  two  types  are  all  valid.  The  basis  functions  over  a  regular  element 
can  be  directly  obtained  from  Eq.(2).  For  an  irregular  element,  Eq.  (7)  can  be  used  to  calculate  the  basis  functions 
with  support  on  it.  Eq.  (18)  defines  the  generalized  Catmull-Clark  basis  functions  over  an  invalid  element,  any  node 
of  which  can  be  an  extraordinary  node.  The  generalized  Catmull-Clark  basis  functions  also  satisfy  partition  of  unity, 
the  proof  of  which  is  the  same  as  that  in  [28]  except  that  generalized  subdivision  matrix  S  is  involved.  However,  we 
do  not  allow  all  its  four  corners  to  be  valence-3,  in  which  case  the  basis  functions  are  linearly  dependent  on  the  invalid 
element  [18].  Moreover,  Eq.  (13)  indicates  refinability  is  also  valid  for  generalized  Catmull-Clark  basis  functions. 
Therefore  we  can  use  them  to  construct  eTHCCS.  With  the  generalized  Catmull-Clark  basis  functions,  preprocessing 
of  input  control  meshes  is  no  longer  required  to  refine  invalid  quadrilaterals,  which  is  a  significant  improvement  on 
efficient  local  refinement,  especially  for  complex  quadrilateral  meshes  with  many  extraordinary  points. 


3.2.  Basis-Function-Insertion  Scheme 

The  original  development  of  THCCS  (or  HB-splines,  THB-splines)  employs  a  basis-function-refinement  scheme, 
which  replaces  a  basis  function  with  its  children  to  enlarge  the  spline  space.  As  a  result,  we  need  to  refine  all  the 
elements  within  the  support  of  to-be-refined  basis  functions.  In  cubic  splines  and  Catmull-Clark  subdivision,  this 
leads  to  the  refinement  of  all  the  two-ring  neighboring  elements,  which  is  not  efficient  to  capture  abrupt  change  in 
geometric  or  solution  features.  Instead,  in  the  following  we  develop  a  new  basis-function-insertion  scheme.  With 
this  scheme,  we  only  need  to  select  the  support  of  inserted  high-level  basis  functions  as  the  to-be-refined  region.  In 
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contrast,  the  to-be-refined  region  in  THB- splines  and  THCCS  is  the  support  of  low-level  basis  functions.  For  example 
we  consider  cubic  splines  in  Fig.  4.  A  basis  function  at  Level  £  has  support  over  the  two-ring  neighboring  Level-/’ 
elements,  as  shown  in  blue  in  Fig.  4(a).  A  Level-(^  +1)  basis  function  has  support  over  the  one-ring  neighboring 
Level-/’  elements  (Fig.  4(b)),  and  a  Level-(/’  +  2)  basis  function  has  support  only  on  one  Level-/1  element  (Fig.  4(c)). 
In  the  following,  we  study  how  to  insert  Level-(/’  +1)  basis  functions,  which  requires  refinement  of  their  one-ring 
neighboring  elements  only. 


- v. 

(c)  Level  £  +  2 


Figure  4:  The  support  of  basis  functions. 


(a) 


(b) 


(c) 


(g) 


(h) 


(i)  0) 


Figure  5:  The  identification  of  to-be-refined  one-ring  neighboring  elements  of  a  regular  node  (a),  a  valence-3  extraordinary  node  (b),  two  nodes  (c), 
and  four  nodes  (d).  Refinement  of  one-ring  neighboring  elements  of  a  regular  node  (f),  a  valence-3  extraordinary  node  (g),  two  nodes  (h),  and  four 
nodes  (i).  (e)  and  (j)  are  equivalent  case  of  the  basis-function-refinement  scheme. 


In  the  basis-function-insertion  scheme,  at  each  refinement  step  we  need  to  identify  and  refine  the  one-ring  neigh¬ 
boring  elements  of  one  or  multiple  nodes.  Local  refinement  can  be  triggered  by  a  user-defined  region  or  simulation 
error.  In  geometric  design,  designers  locally  refine  regions  of  interest  to  add  more  features.  In  analysis,  elements  with 
large  error  need  to  be  refined  to  improve  numerical  performance.  Assume  that  we  use  simulation  error  to  identify 
to-be-refined  elements.  We  group  all  the  one-ring  neighboring  elements  for  each  Level-/’  node  and  compare  their 
error  with  a  given  threshold.  Then  a  set  of  elements  is  identified  as  to-be-refined.  In  the  control  mesh  we  have  both 
regular  nodes  and  extraordinary  nodes.  Let  us  first  study  the  refinement  of  one-ring  neighboring  elements  of  a  regular 
node.  Assume  we  have  a  local  mesh  at  Level  £  {£  >  0)  shown  in  Fig.  5(a),  where  a  regular  node  is  marked  with  a  red 
dot  and  its  four  one-ring  neighboring  elements  are  to  be  refined,  as  marked  in  blue.  After  refinement,  we  obtain  16 
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Level-(£  +1)  elements  as  shown  in  Fig.  5(f).  Next  we  define  the  basis-function-insertion  rule  to  insert  a  Level-(^  +1) 
basis  function. 

Basis-Function-Insertion  Rule.  During  local  refinement  around  a  node  at  Level  £,  a  Level-(£  +  1 )  basis  function  has 
to  be  inserted  associated  with  this  node  if  all  its  two-ring  Level-(£  +  1)  neighboring  elements  are  generated. 

The  basis-function-insertion  rule  is  straightforward  since  a  (generalized)  Catmull-Clark  basis  function  has  support 
over  its  two-ring  neighboring  elements.  Back  to  Fig.  5(f),  the  newly  generated  16  Level- (f  +1)  elements  are  the 
two-ring  neighborhood  of  the  green  dot.  Therefore  according  to  the  basis-function-insertion  rule,  there  exists  one 
basis  function  associated  with  this  node.  In  this  manner,  the  refinement  of  the  one-ring  neighboring  elements  of  a 
regular  node  at  Level  £  leads  to  the  insertion  of  a  Level-(^  +1)  basis  function. 

The  same  idea  can  be  used  to  handle  extraordinary  nodes.  In  Fig.  5(b),  the  elements  in  blue  are  the  one-ring 
neighboring  elements  of  a  valence- 3  extraordinary  node  (red  dot).  They  are  refined  and  12  Level- {£  +1)  elements 
are  generated,  as  shown  in  Fig.  5(g).  According  to  the  basis-function-insertion  rule,  the  Level-(^  +1)  basis  function 
associated  with  the  green  dot  is  inserted.  In  general,  for  the  extraordinary  node  of  any  valence  N ,  by  refining  its 
one-ring  neighboring  elements  we  obtain  4 N  Level-(^  +1)  elements  together  with  one  inserted  Level-(^  +1)  basis 
function. 

The  number  of  inserted  high-level  basis  functions  varies  when  the  one-ring  neighboring  elements  of  multiple 
nodes  are  refined.  For  instance,  Fig.  5(c)  shows  the  one-ring  neighboring  elements  (in  blue)  of  two  nodes  of  an  edge. 
The  refinement  of  these  elements  is  shown  in  Fig.  5(h),  where  we  observe  that  3  basis  functions  associated  with  the 
green  dots  have  to  be  inserted  according  to  the  basis-function-insertion  rule.  Moreover,  in  the  case  of  four  corner 
nodes  of  an  element  as  shown  in  Fig.  5(d),  9  Level-(£  +1)  basis  functions  are  inserted;  see  Fig.  5(i).  When  we 
have  the  case  in  Fig.  5(e),  the  one-ring  neighboring  elements  of  those  red  dots  are  actually  the  two-ring  neighboring 
elements  of  the  valence-3  extraordinary  node.  Thus  the  refinement  leads  to  the  equivalent  case  of  the  basis-function- 
refinement  scheme.  The  basis  function  of  the  extraordinary  node  is  replaced  by  its  children  associated  with  green  dots 
in  Fig.  5(f).  In  general,  the  basis-function-insertion  scheme  yields  less  refinement  than  the  basis-function-refinement 
scheme.  Practical  cases  can  be  more  complicated,  but  we  can  always  follow  the  basis-function-insertion  rule  to 
determine  where  the  basis  functions  need  to  be  inserted.  The  insertion  of  the  basis  function  enlarges  the  spline  space, 
leading  to  the  nested  property,  which  will  be  proved  later  in  Appendix  B.  The  basis-function-refinement  scheme  can 
also  be  directly  applied  to  truncated  hierarchical  B- splines. 


3.3.  Truncation 


Figure  6:  Five  examples  of  Level-/7  to-be-truncated  basis  functions  (green  circles). 


Inserting  high-level  basis  functions  destroys  partition  of  unity  and  changes  the  geometry.  To  resolve  this  issue, 
we  apply  a  truncation  mechanism  to  the  neighboring  low-level  basis  functions.  A  basis  function  needs  to  be  truncated 
if  any  of  its  children  is  added  in  the  spline  space.  Let  &  be  the  set  of  Level-/’  basis  functions.  The  inserted  basis 
functions  are  at  Level  £  +  1  and  they  are  added  into  the  spline  space,  denoted  as  H3M.  Then  some  Level-/’  basis 
functions  (&t)  need  to  be  truncated  if  any  of  their  children  is  an  inserted  Level-(/’  +1)  basis  function,  and  we  identify 
these  basis  functions  as 


=  {B\  e  |  chd B\  n  ±  0}. 


(20) 
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For  the  five  cases  studied  in  Fig.  5,  Fig.  6  illustrates  how  Eq.  (20)  is  used  to  identify  to-be-truncated  basis  functions. 
After  refining  the  blue  elements,  we  need  to  truncate  basis  functions  associated  with  the  corner  nodes  of  all  the  refined 
elements,  as  marked  with  green  circles.  According  to  the  refinability  property,  a  Level-^  basis  function  Bl.  can  be 
expressed  as  a  linear  combination  of  its  children,  and  we  have 


B\=  2  C'jB1;1,  (21) 

B^'  echdflf 

where  comes  from  the  general  subdivision  matrix  S  or  A.  The  truncation  is  performed  for  basis  functions  Bf  e  Bft 
by  removing  the  children  contained  in  B(  from  the  summation  in  Eq.  (21),  that  is, 


trun£f=  Yj  c‘JbY-  (22) 

B J+1  echdfif  and  Bf+]  £Bi+x 

Note  that  in  Fig.  6(e),  all  the  children  of  the  basis  function  associated  with  the  blue  dot  are  inserted  Level-(£  +1)  basis 
functions;  see  Fig.  5(j).  Therefore  this  basis  function  becomes  a  zero  function  after  truncation.  In  this  case  we  call 
such  basis  functions  passive  and  they  no  longer  exist  in  the  spline  space.  The  other  non-zero  basis  functions  are  called 
active.  The  passive  basis  functions  are  actually  the  refined  basis  functions. 

Truncating  basis  functions  can  help  us  remove  the  repeated  contribution  of  high-level  basis  functions.  Consider 
a  Level-^  element  £lk  and  let  Ilk  be  the  index  set  of  all  Level-^  basis  functions  Bl{  with  support  on  it.  Suppose  before 
inserting  any  new  Level-(^  +1)  basis  functions,  Bl  satisfy  partition  of  unity  and  we  have  Z; ei{  B£  =  1.  According  to 
refinability,  we  further  obtain  Tuiei*  B f  =  ZfejJ  Yjjece  cijBfj+1  =  k  where  C\  is  the  index  set  of  children  basis  functions 
of  Bfr  Now  we  insert  a  Level-(^  +1)  basis  function  ( B ^+1)  by  refinement  and  B ^,+1  has  support  over  Q?k:  The  summation 
of  basis  functions  over  Qfk  becomes  ZyeC  cijBi]+X  +  B^1  =  1  +  Z^,+1.  This  can  be  rewritten  as 


M 

j 


+  B^1 


z  z 

ieIl 


+  bS+1 


Z  Z  cijB';'  +  2B';\ 

ieIl  JeCt’frf 


(23) 


where  Z  iei*  cijf  B^?1  =  B^]  Z/e/f  Gy  =  B^1  because  Z/G/f  cif  =  1  holds  for  Catmull-Clark  subdivision  and  knot 
insertion  algorithm  [28].  From  Eq.  (23)  we  can  observe  that  the  summation  of  basis  functions  over  £lk  counts  the 
inserted  basis  function  B ^,+1  twice.  By  removing  CiyB ^+1  from  the  refinability  relationship,  we  can  achieve  partition  of 
unity.  Furthermore,  the  geometry  is  preserved  during  local  refinement,  which  will  be  proved  in  Appendix  B. 


3.4.  Construction  ofeTHCCS 

Similar  to  the  construction  of  THCCS  [28],  we  follow  three  steps  to  construct  eTHCCS:  identification  of  to-be- 
refined  elements  and  to-be-truncated  basis  functions  (Step  1),  refinement  of  identified  elements  and  truncation  of 
identified  basis  functions  (Step  2),  and  collection  of  all  the  hierarchical  basis  functions  and  elements  (Step  3). 

We  start  with  an  initial  control  mesh  that  can  be  any  quadrilateral  mesh,  except  the  one  with  all  its  vertices  that 
are  valence-32.  Given  the  input  control  mesh  At0,  the  initial  set  of  eTHCCS  basis  functions  is  defined  by  Level-0 
basis  functions  B°  associated  with  At0,  and  the  entire  domain  is  defined  as  suppiB0.  The  initial  elements  of  all  the 
quadrilaterals  in  At0  are  denoted  as  8°.  eTHCCS  is  recursively  constructed  up  to  a  desired  maximum  level,  ^max 
(^max  >  0),  which  allows  us  to  study  two  consecutive  levels  (Level  €  and  Level  €  +  1)  at  one  time.  Suppose  we  have 
constructed  Level-^  elements  and  basis  functions,  and  now  we  want  to  construct  Level  €  +  1.  Let  &  be  the  set  of 
Level-^  basis  functions  and  8£  be  the  set  of  Levels  elements.  8l  defines  the  sub-domain  (Q^)  at  Level  l.  The  eTHCCS 
basis  functions  of  l  levels  are  collected  in  the  set  ^fTHCCS,  whereas  the  eTHCCS  elements  are  in  £fXHCCS- 

Identification  (Step  1).  As  discussed  in  Section  3.2  we  use  simulation  error  to  identify  to-be-refined  elements. 
Then  a  set  of  elements  is  identified  as  to-be-refined  if  their  error  is  larger  than  a  given  threshold,  denoted  as  8lr. 
Besides,  all  the  basis  functions  associated  with  the  corner  nodes  of  these  elements  are  identified  as  to-be-truncated 
basis  functions  ( &t);  see  Fig.  6(a-e). 


2  A  mesh  with  all  valence-3  vertices  produces  linearly  dependent  blending  functions  [18]  and  thus  cannot  be  used  in  analysis. 
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Refinement  and  Truncation  (Step  2).  Refinement  of  elements  in  can  be  easily  obtained  by  quadtree  subdi¬ 
vision.  After  refinement,  elements  in  are  defined  as  passive  whereas  the  newly  generated  Level-(^  +1)  elements 
(SM)  are  defined  as  active.  Only  active  elements  will  be  used  in  eTHCCS  construction.  All  the  Level-(^  +1)  elements 
define  the  sub-domain  at  Level  l  +  1,  denoted  as  Qf+l.  It  is  obvious  that  c  Of.  We  check  which  node  has  all  its 
two-ring  neighboring  elements  generated.  According  to  the  basis-function-insertion  rule,  Level-(^  +1)  basis  functions 
(tBM)  are  inserted  for  such  nodes  and  they  will  be  added  in  the  eTHCCS  space.  The  corresponding  Level-(^  +1)  con¬ 
trol  points  are  calculated  by  the  Catmull-Clark  subdivision  rule.  Next,  the  Level-^  basis  functions  in  &t  are  truncated 
according  to  Eq.  (22).  With  all  its  children  inserted,  a  basis  function  becomes  zero  and  it  is  passive,  so  it  no  longer 
exists  in  eTHCCS. 

Collection  (Step  3).  On  one  hand,  by  refinement  elements  in  become  passive  and  basis  functions  in  are 
truncated,  some  of  which  are  even  passive  ones.  Therefore,  we  remove  elements  from  £fTHCCS  and  update  the 
basis  functions  of  &t  in  ^XHCCS-  On  the  other  hand,  we  obtain  new  active  Level-(^  +1)  elements  and  basis  functions, 
which  are  used  to  construct  eTHCCS  of  (£  +  1)  levels.  We  have 


and 


pM 

°eTHCCS 


<Rp+ 1 

^eTHCCS 


_  pt  1  I  pt+ 1 

~~  °eTHCCS  u  ° 

(24) 

-  U  &+l 

~  ^eTHCCS  u 

(25) 

We  can  recursively  perform  Step  1  to  Step  3  until  the  maximum  level  £max  is  reached.  The  construction  enlarges 
the  spline  space  of  eTHCCS  with  nested  sub-domains  as  the  level  increases,  that  is, 

n°  2  fl1  2  •  •  •  3  Q4-”  (26) 

and 

span^eTHCCS  c  span^eTHCCS  c  •  •  •  c  span^SeTHCCS.  (27) 

We  will  prove  this  property  in  Appendix  B. 

Remark  3.2.  Given  a  local  region  at  Level  A  in  this  paper  we  focus  on  inserting  Level-(^  +1)  basis  functions  such 
that  we  can  recursively  construct  eTHCCS  level  by  level.  However,  inserting  basis  functions  at  higher  levels  is  also 
supported  in  our  algorithm.  For  instance,  inserting  a  Level-(^  +  2)  basis  function  requires  refinement  of  a  Level-^ 
element  into  16  Level-(^  +  2)  elements,  and  the  truncation  of  a  Level-^  basis  function  needs  its  Level-(^  +  2)  children 
basis  functions.  The  construction  procedure  is  the  same  as  discussed  in  Sections  3.2  and  3.3. 


4.  Examples  and  Discussion 

In  this  section,  we  first  study  the  analysis-suitability  of  generalized  Catmull-Clark  basis  functions  via  three  patch 
tests.  Then  we  solve  a  benchmark  problem:  the  Laplace  equation  on  the  L-shaped  domain,  where  both  THCCS  [28] 
and  eTHCCS  basis  functions  are  used  to  study  the  convergence  behavior.  In  the  end,  we  solve  the  Laplace  equation 
over  four  complex  models.  Table  1  summarizes  the  statistics  of  the  tested  models. 

4.1.  Patch  Tests 

For  patch  tests,  we  solve  a  2D  linear  elasticity  problem  by  applying  uniform  tension  to  a  square  (Young’s  mod¬ 
ulus  E  -  1,  Poisson’s  ratio  /i  =  0.3).  We  use  three  input  irregular  quadrilateral  meshes  with  different  numbers  of 
extraordinary  nodes,  as  shown  in  Figs.  7(a)  ~  9(a).  The  meshes  in  8(a)  and  9(a)  have  elements  with  more  than  one 
extraordinary  node.  In  particular,  all  the  nodes  of  the  central  element  are  valence-3  in  9(a).  Recall  that  according  to 
Peters  [18],  subdivision  functions  are  not  linearly  independent  on  such  elements,  which,  however,  does  not  violate  the 
global  linear  independence  condition  unless  all  the  nodes  in  the  mesh  are  valence-3.  Our  application  only  requires 
global  linear  independence,  and  generalized  Catmull-Clark  basis  functions  are  used  here.  To  accurately  integrate  an 
element  with  extraordinary  nodes,  we  subdivide  the  elements  into  a  sequence  of  sub-elements,  as  indicated  in  Fig. 
1(d),  and  in  each  sub-element,  we  adopt  4x4  Gaussian  integration.  Two  strain  components  in  each  patch  test  are 


11 


Table  1 :  Statistics  of  all  the  tested  models 


Models 

#  Input 
Nodes 

#  Input 
Elements 

#  Invalid 

Elements 

£ 

#  Levels 

THCCS 

#  Levels 

eTHCCS 

DOF- 

THCCS 

DOF- 

eTHCCS 

DOF 

Ratio 

L-Regular 

153 

128 

0 

2E-4 

5 

7 

773 

391 

50.6% 

L-Irregular 

221 

192 

0 

2.5E-3 

4 

6 

807 

549 

68.0% 

Genus- 3 

3,068 

3,072 

0 

0.1 

5 

5 

4,016 

3,427 

85.3% 

Bunny 

3,023 

3,021 

4 

0.1 

5 

5 

4,682 

3,482 

74.4% 

Venus 

1,559 

1,552 

899 

0.5 

4 

5 

6,646 

1,822 

27.4% 

Head 

2,909 

2,907 

1,572 

0.5 

5 

4 

11,318 

3,207 

28.3% 

Note:  #  stands  for  number  and  s  is  the  given  threshold.  DOF  Ratio  =  (eTHCCS  DOF)/(THCCS  DOF). 


shown  in  Figs.  7(b,  c)  ~  9(b,  c),  where  the  black  curves  are  the  isoparametric  lines  projected  onto  the  physical  do¬ 
main.  We  calculate  the  error  in  the  L 2  norm  and  Hl  norm,  and  display  them  with  respect  to  subdivision  levels  used  for 
Gauss  integration  of  irregular  elements;  see  Figs.  7(d)  ~  9(d).  The  error  decreases  as  the  subdivision  level  increases. 
However,  it  remains  of  the  same  order  when  we  use  more  than  15  levels.  As  discussed  in  [20,  17,  28],  Catmull-Clark 
basis  functions  have  limitations  in  analysis.  Directly  applying  Gaussian  quadrature  over  elements  with  extraordinary 
nodes  introduces  numerical  error,  because  Catmull-Clark  basis  functions  are  infinite  piecewise  polynomials.  Further 
study  is  needed  to  develop  efficient  and  accurate  quadrature  schemes. 


(a)  (b)  (c)  (d) 


Figure  7:  Patch  test  1  on  an  irregular  quadrilateral  mesh  with  2  extraordinary  nodes  using  generalized  Catmull-Clark  basis  functions,  (a)  The  input 
mesh  and  boundary  conditions;  (b,  c)  two  strain  components  in  X  -  X  and  Y  -  Y  directions,  respectively;  and  (d)  error  with  respect  to  subdivision 
levels  used  for  Gauss  integration. 


(a) 


Figure  8:  Patch  test  2  on  an  irregular  quadrilateral  mesh  with  6  extraordinary  nodes  using  generalized  Catmull-Clark  basis  functions,  (a)  The  input 
mesh  and  boundary  conditions;  (b,  c)  two  strain  components  in  X  -  X  and  Y  -  Y  directions,  respectively;  and  (d)  error  with  respect  to  subdivision 
levels  used  for  Gauss  integration. 


4.2.  L- shaped  Problem 

Fig.  10(a)  shows  the  problem  setting  of  the  Laplace  equation  Au  =  0  over  the  L-shaped  domain  [-1, 1]2\[0,  l]2, 
where  the  Dirichlet  boundary  conditions  (T D)  are  strongly  imposed.  The  analytical  solution  is  available  in  polar 


12 


Figure  9:  Patch  test  3  on  an  irregular  quadrilateral  mesh  with  8  extraordinary  nodes  using  generalized  Catmull-Clark  basis  functions,  (a)  The  input 
mesh  and  boundary  conditions;  (b,  c)  two  strain  components  in  X  -  X  and  Y  -  Y  directions,  respectively;  and  (d)  error  with  respect  to  subdivision 
levels  used  for  Gauss  integration. 


Figure  10:  Laplace  equation  on  the  L-shaped  domain,  (a)  Geometry  and  problem  settings;  (b)  input  regular  control  mesh;  and  (c)  input  irregular 
control  mesh. 


coordinates  (r,  0), 


u(r ,  0)  =  r2/3  sin(2#/3  -  n/3),  where  r  >  0  and  7r/2  <  6  <  2n.  (28) 

We  use  two  input  control  meshes:  a  regular  mesh  and  an  irregular  mesh,  shown  in  Fig.  10(b,  c),  respectively.  For  each 
mesh,  three  refinement  schemes  are  studied:  uniform  refinement,  THCCS  [28]  and  eTHCCS.  The  uniform  refinement 
simply  subdivides  all  the  elements  at  each  refinement  step.  The  error  is  assessed  in  the  L 2  norm  and  Hl  semi-norm 
for  each  element,  as  well  as  the  entire  domain.  In  THCCS  or  eTHCCS,  to-be-refined  elements  are  identified  in  terms 
of  two-ring  or  one-ring  neighborhood  of  a  node.  Therefore,  we  convert  the  element-wise  error  to  the  node-wise 
error  by  summing  the  error  of  two-ring  (in  THCCS)  or  one-ring  (in  eTHCCS)  neighboring  elements  of  the  node.  At 
each  refinement  step,  we  refine  elements  with  the  node- wise  error  larger  than  rj  •  emax,  where  rj  (0%  <  rj  <  100%) 
refers  to  a  prescribed  percentage  and  emax  is  the  maximum  node- wise  error.  In  this  problem,  we  set  rj  =  100%  to 
refine  the  elements  with  maximum  node-wise  error.  To  improve  the  accuracy  of  numerical  integration  surrounding 
an  extraordinary  node,  if  any  element  within  its  two-ring  neighborhood  is  to  be  refined,  we  refine  all  its  two-ring 
neighboring  elements.  The  adaptive  analysis  terminates  when  the  L 2  error  (or  Hl  semi-norm)  over  the  entire  domain 
is  smaller  than  a  given  threshold  s.  Fig.  11  (a,  c)  shows  the  distribution  of  element-wise  L 2  error  using  THCCS  at 
the  final  step,  and  Fig.  1 1  (b,  d)  shows  this  result  using  eTHCCS.  We  observe  that  the  refinement  on  both  regular  and 
irregular  meshes  is  more  localized  at  the  sharp  corner,  where  there  is  a  singularity  in  the  solution.  In  Fig.  12,  the  L 2 
error  and  Hl  semi-norm  are  plotted  with  respect  to  degrees  of  freedom  (DOF).  eTHCCS  achieves  the  same  accuracy 
with  only  50.6%  DOF  of  THCCS  in  the  regular  mesh  and  68%  of  THCCS  in  the  irregular  mesh. 


13 


(a)  (b)  (c)  (d) 


Figure  11:  L2  error  distribution  of  numerical  solutions  on  the  L-shaped  domain  at  the  final  refinement  step:  (a)  the  regular  mesh  with  THCCS;  (b) 
the  regular  mesh  with  eTHCCS;  (c)  the  irregular  mesh  with  THCCS;  and  (d)  the  irregular  mesh  with  eTHCCS. 


4.3.  Complex  Models 

We  also  solve  the  Laplace  equation  on  four  complex  models:  the  genus- 3  model  (Fig.  13),  the  bunny  model  (Fig. 
14),  the  Venus  model  (Fig.  15)  and  the  head  model  (Fig.  16).  The  input  quadrilateral  meshes  have  no  elements  with 
nodes  that  are  all  valence-3.  However,  the  bunny  model,  the  Venus  model  and  the  head  model  have  invalid  elements, 
where  generalized  Catmull-Clark  basis  functions  are  used.  To  create  a  solution  field  with  abrupt  local  change,  we 
strongly  prescribe  the  solution  to  have  a  certain  value,  say  100,  on  some  elements,  and  then  set  a  very  different  value, 
say  1,  over  some  elements  nearby.  And  we  study  the  performance  of  THCCS  and  eTHCCS.  Figs.  13(a)  -  16(a) 
show  the  input  meshes  with  boundary  conditions,  where  the  elements  marked  in  red  are  set  at  100  and  the  elements 
in  blue  are  set  at  1.  Due  to  the  lack  of  analytical  solutions,  the  L 2  error  for  each  element  is  approximated  by  using  the 
so-called  bubble  function  [25].  Then  following  the  same  procedure  as  in  solving  the  L-shaped  problem,  we  perform 
adaptive  analysis  until  the  relative  maximum  error  is  smaller  than  a  given  threshold  s.  The  error  is  defined  as  en /e®, 
where  en  is  the  maximum  element-wise  error  after  n  refinement  steps  and  e°  is  the  maximum  element-wise  error  of 
THCCS  at  the  initial  step.  We  set  s  =  0.1  for  the  genus-3  and  bunny  models,  and  s  =  0.5  for  the  Venus  and  head 
models.  We  also  set  rj  =  70%  to  refine  more  elements  than  for  the  L-shaped  problem  at  each  refinement  step. 

At  the  final  step,  the  solutions  using  THCCS  are  shown  in  Figs.  13(b)  -  16(b)  whereas  the  solutions  using 
eTHCCS  are  shown  in  Figs.  13(c)  -  16(c).  From  Table  1,  eTHCCS  is  the  most  efficient  method,  especially  in  the 
Venus  and  head  models,  where  at  the  final  step,  the  degrees  of  freedom  are  only  27.4%  and  28.3%  of  those  using 


(a) 


(b) 


Figure  12:  Convergence  curves  with  respect  to  L 2  error  (a)  and  H]  error  (b).  The  legend  “UR”  represents  uniform  refinement. 
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Figure  13:  Solving  Laplace  equation  over  a  genus-3  model,  (a)  Input  quadrilateral  mesh  with  boundary  conditions;  (b)  numerical  solution  using 
THCCS;  (c)  numerical  solution  using  eTHCCS;  (d,  e)  zoom-in  picture  of  the  window  in  (b,  c)  respectively. 


Figure  14:  Solving  Laplace  equation  over  a  bunny  model,  (a)  Input  quadrilateral  mesh  with  boundary  conditions;  (b)  numerical  solution  using 
THCCS;  (c)  numerical  solution  using  eTHCCS;  (d,  e)  zoom-in  picture  of  the  window  in  (b,  c)  respectively. 
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THCCS.  The  significant  improvement  mainly  benefits  from  the  generalized  Catmull-Clark  basis  functions.  As  the 
input  quadrilateral  meshes  of  the  Venus  model  and  the  head  model  have  a  large  number  of  invalid  elements,  THCCS 
needs  to  refine  all  such  elements,  resulting  in  an  almost  uniform  refinement;  see  Figs.  15(b)  and  16(b).  In  contrast, 
generalized  Catmull-Clark  basis  functions  can  be  directly  applied  to  the  input  meshes  and  the  refinement  is  only 
performed  for  elements  with  large  error,  as  shown  in  Figs.  15(c)  and  16(c).  On  the  other  hand,  the  genus-3  model 
and  the  bunny  model  only  have  a  few  invalid  elements.  The  improvement  of  efficiency  mainly  benefits  from  the 
basis-function-insertion  scheme.  eTHCCS  only  uses  85.3%  DOF  of  THCCS  in  the  genus-3  model  and  74.4%  DOF  of 
THCCS  in  the  bunny  model.  The  basis-function-insertion  scheme  works  well  especially  when  the  solution  field  has 
significant  local  features. 


5.  Conclusion  and  Future  Work 

In  conclusion,  we  develop  eTHCCS  with  generalized  Catmull-Clark  basis  functions  and  a  new  basis-function- 
insertion  scheme,  aiming  to  improve  the  efficiency  of  local  refinement.  The  generalized  Catmull-Clark  basis  functions 
directly  work  on  invalid  elements  with  more  than  one  extraordinary  node,  providing  a  basis  for  isogeometric  analysis 
of  arbitrary  quadrilateral  meshes.  The  basis-function-insertion  scheme  releases  the  selection  restriction  of  the  to-be- 
refined  region  in  THB- splines  and  THCCS.  In  cubic  splines,  it  allows  refinement  of  a  single  element  at  each  refinement 
step,  while  THB-splines  or  THCCS  need  to  refine  at  least  two-ring  neighboring  elements.  In  practice,  as  we  construct 


Figure  15:  Solving  Laplace  equation  over  a  Venus  model,  (a)  Input  quadrilateral  mesh  with  boundary  conditions;  (b)  numerical  solution  using 
THCCS;  (c)  numerical  solution  using  eTHCCS;  (d,  e)  zoom-in  picture  of  the  window  in  (b,  c)  respectively. 
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Figure  16:  Solving  Laplace  equation  over  a  head  model,  (a)  Input  quadrilateral  mesh  with  boundary  conditions;  (b)  numerical  solution  using 
THCCS;  (c)  numerical  solution  using  eTHCCS;  (d,  e)  zoom-in  picture  of  the  window  in  (b,  c)  respectively. 


eTHCCS  level  by  level,  we  refine  one-ring  neighboring  elements.  The  basis-function-insertion  scheme  can  also  be 
applied  to  truncated  hierarchical  B-splines.  Like  THCCS,  eTHCCS  preserves  geometry  and  produces  nested  spline 
spaces.  eTHCCS  basis  functions  also  satisfy  partition  of  unity,  convex  hull  and  global  linear  independence.  We  use 
five  numerical  models  to  test  the  proposed  method.  From  the  results,  we  can  observe  that  eTHCCS  achieves  the  same 
accuracy  with  fewer  degrees  of  freedom  than  THCCS.  In  the  future,  it  may  be  promising  to  employ  the  proposed 
method  on  hierarchical  T- splines,  where  the  patch  test  can  be  successfully  passed  without  introducing  huge  numbers 
of  quadrature  points. 
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Appendix  A:  Generalized  Subdivision  Matrix  and  Selection  Matrices 

For  an  element  with  more  than  one  extraordinary  nodes,  we  develop  the  generalized  subdivision  matrix  to  perform 
local  Catmull-Clark  subdivision.  The  generalized  subdivision  matrix  S j3Xi3  is  used  to  calculate  the  new  y'3  vertices 
at  Level  £  +  1  (Fig.  2(b))  from  the  neighboring  z'3  vertices  at  Level  l  (Fig.  2(a)).  S  is  constructed  directly  by  the 
Catmull-Clark  subdivision  rule.  Following  the  indices  locally  labeled  as  in  Fig.  2,  we  have 
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where  aNt  =  l  -  - h,bNk  =  ^,cNt  =  jL  (k  =  0, 1,2,3),  d  =  §,  e  =  and  /  =  The  generalized  subdivision 

’ k  yk  yk 

matrix  depends  on  the  valence  number  of  four  corners  of  the  element. 


Appendix  B:  Geometry  Preservation  and  Nested  Property 

In  Section  3,  we  claim  that  eTHCCS  can  preserve  the  geometry  and  possesses  nested  property.  In  the  following, 
we  mathematically  prove  these  two  properties. 

Proposition  1.  The  geometry  is  preserved  during  the  construction  of  eTHCCS  from  Level  £  (£  >0)  to  Level  7  +  1. 

Proof  We  prove  this  proposition  by  constructing  Level  7  +  1  from  Level  £  (7  >  0),  showing  that  the  geometry  is  the 
same  before  and  after  refinement.  After  refinement,  active  Level-7  elements  remain  the  same,  leading  to  the  same 
geometry  as  before.  Therefore,  we  only  need  to  focus  on  the  Level-(7  +1)  elements.  Let  Qfk+l  be  a  Level-(7  +  1) 
element  obtained  by  refining  a  Level-7  element  Thus  we  have  Qfk+l  c  Q?  After  refinement,  the  portion  of  the 
geometry  (*S|q£+i)  corresponding  to  Q£+1  is  calculated  as 

SIqm  =  L  C'B?'  +  L  7trull/f  +  E  74  (30) 

iel%+l  jeT{  jeI(a\T{ 

where  7^,  I£a+l  denote  the  index  set  of  active  Levels  and  Level-(7  +1)  non-truncated  basis  functions,  Tl  is  the  index 
set  of  Levels  truncated  basis  functions,  and  Pj,  P*+l  are  Level-7  and  Level-(7  +1)  control  points,  respectively.  Eq. 
(30)  consists  of  a  summation  of  three  terms  because  active  Level-(7  +1)  basis  functions  (7^+1),  Level-7  truncated  basis 
functions  (trun^)  and  other  active  Level-7  basis  functions  (B\)  may  all  have  support  on  Cl[+l .  According  to  Eq.  (22), 
trun Bi  can  be  expressed  as 

trunfK  =  Yj  (31) 

ieIM  \Ta+] 

where  IM  represents  an  index  set  of  basis  functions  associated  with  the  Catmull-Clark  mesh  obtained  by  7  +  1 
subdivisions.  Note  that  an  active  non-truncated  Level-7  basis  function  Be.  ( j  e  I{a\T{)  does  not  have  any  active 
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children  at  Level  £+\  (otherwise  it  contradicts  the  definition  of  truncation).  According  to  Eq.  (21),  we  have 


bj=  Z  c^+1- 

ieIM  \/f 

Substituting  Eqs.  (31)  and  (32)  into  Eq.  (30),  we  have 
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(33) 


Note  that  Z7e/f  ci^j  ~  ^f+1  (*  G  /+1\/^+1)  directly  comes  from  the  Catmull-Clark  subdivision  rule.  Recall  that  the 
limit  Catmull-Clark  subdivision  surface  can  be  equivalently  calculated  from  any  control  mesh  in  the  global  refinement 
sequence.  Thus,  Eq.  (33)  means  that  the  limit  surface  is  calculated  by  a  Level-(^  +1)  Catmull-Clark  control  mesh. 
Any  other  portions  of  the  geometry  can  be  handled  in  the  same  manner.  Therefore,  the  limit  surface  does  not  change 
during  the  construction  of  eTHCCS. 

□ 


The  proof  of  Proposition  1  is  similar  to  the  proof  of  geometry  preservation  in  [28].  In  Eq.  (31),  Cji  can  all  be 
zero  and  Bl.  does  not  contribute  to  the  calculation  of  <S|q^+i  .  Actually,  this  case  implies  that  the  children  of  Bl.  are 
all  active  basis  functions  at  Level  £  +  1  and  Bl.  becomes  passive  at  Level  £,  which  is  exactly  the  same  case  as  for 
THCCS  construction  [28].  The  eTHCCS  space  also  has  nested  property  as  described  in  Eq.  (26)  and  Eq.  (27).  As 
the  hierarchical  level  increases,  the  sub-domain  decreases  and  the  eTHCCS  is  enlarged.  Since  eTHCCS  is  recursively 
constructed,  without  loss  of  generality,  we  only  need  to  prove  that  the  nested  property  holds  for  two  consecutive  levels, 
Level  £  {£  >  0)  and  Level  £  +  1 . 


Proposition  2.  Given  an  eTHCCS  with  levels  up  to  £  (£  >  0),  Level  £+  l  is  constructed  using  the  basis -function- 
insertion  scheme.  The  eTHCCS  space  is  enlarged,  that  is, 


spantB y 


eTHCCS 


c  spantBfe+l 


eTHCCS ’ 


(34) 


where  &eTHCCS  and  &eTHCCS  contain  the  eTHCCS  basis  functions  of  £  levels  and  £  +  1  levels,  respectively. 

Proof.  To  prove  Eq.  (34),  we  only  need  to  prove  each  basis  function  in  ^fTHCCS  can  be  represented  by  a  linear 

During  eTHCCS  construction,  the  to-be-truncated  Level-^  basis  functions 

The  other  basis  functions  in 


combination  of  basis  functions  in  ^thccs  • 

in  ^fTHCCS  are  used  to  construct  the  Level-^  truncated  basis  functions  in  ^xhccs- 

^eTHCCS  remain  ^e  same  in  thccs •  Therefore,  we  only  need  to  check  those  to-be-truncated  basis  functions.  Before 
truncation,  a  to-be-truncated  basis  function  B f  e  ^hccs  ' 


can  be  expressed  by  a  linear  combination  of  its  children, 


^  =  ZcA+l  =  ZcA+l  +  Z 

jeCf  jel€+l  jeC.\If+] 


(35) 


where  C\  is  the  index  set  of  the  children  of  B\,  and  r  is  the  index  set  of  newly  inserted  Level-(£  +1)  basis  functions. 
Note  that  Y*iec€\iM  c0'^+1  *s  actually  the  truncated  basis  function  with  respect  to  Bfr  and  we  have 

B\  =  2  cuBep  +  trunflf,  (36) 

jzTa+l 
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where  B ^+1  €  ^xhccs  f°r  J  e  ^+1»  and  trun^  e  ^^hccs*  Therefore,  any  to-be-truncated  basis  function  can  be 
expressed  by  a  linear  combination  of  basis  functions  in  ^thccs*  To  enc*’  we  Prove  that  any  basis  function  in 
traces  can  he  represented  by  a  linear  combination  of  those  in  ^^hccs  •  Therefore  Eq.  (34)  holds  and  Proposition  2 
holds. 

□ 
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